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To investigate the steady thermal hydraulic characteristics of U-tube steam generator (SG), a 1D simulation 
code based on the four-equation drift flux model is developed. The U-tube channels presumably consist mainly 
of the primary channel, secondary channel, and tube wall. In the sub-cooling regions of the primary and sec- 
ondary channels, flow is simulated using the single-phase flow model, whereas that in the boiling regions of the 
secondary channels is simulated using the four-equation drift flux model. The first-order equations of upwind 
difference are derived based on the staggered grid. Steady-state thermal hydraulic parameters are obtained with 
a cross-iteration scheme of heat balance and natural circulation requirement. The developed code is applied to 
analyze the SG behavior of the Qinshan I Nuclear Power Plant under 100%, 75%, 50%, 30%, and 15% power 
conditions. Analysis results are then compared with the simulation results obtained using RELAPS. 


Keywords: U-tube steam generator, Thermal hydraulic characteristic, Steady simulation, Four-equation drift flux model 


DOI: 10.13538/j.1001-8042/nst.25.050601 


I. INTRODUCTION 


The U-tube steam generator (SG) is a heat exchanger that 
connects the primary and secondary coolant loops in a nu- 
clear power plant (NPP). According to worldwide statistics, 
operational accidents of SG-related pressurized water reac- 
tor (PWR) account for a large proportion of all PWR acci- 
dents [1]. Approximately 1/4 of the unplanned outage cases 
in PWR NPPs are caused by SG failure. Given that the 
flow and heat transfer in the primary and secondary loops are 
closely connected to the safety and stable operation of the 
SG, it is of importance to understand their thermal hydraulic 
characteristics. 

The U-tube SG can presumably be a nonlinear, complex 
system with many flowing parameters. Studies on the thermal 
hydraulic behavior of SGs have achieved greatly. The follow- 
ings are examples of SG simulation codes. The THEDA2 
code developed in the U.S. uses 3D conservation equations 
of mass, momentum, and energy for a homogeneous equilib- 
rium mixture (HEM) model [2]. The ATHOS code applies ei- 
ther the three-equation HEM model or the four-equation drift 
flux model with options for 1D, 2D, and 3D analyses [3].The 
THIRST code developed by AECL for 700 MWe SG, is a 
1D thermal hydraulic [4]. Several thermal hydraulic codes 
for NPP SGs were developed in China. Based on the 1D 
separated fluid model, the SGTH-2 performs steady analy- 
sis of the U-tube SG [5].The MOFS is based on the 1D HEM 
model [6], while the SG code for high-temperature gas cool- 
ing reactors follows the 2D HEM model [7]. 

Most of the existing thermal hydraulic codes for NPP SGs 
utilize the classical HEM model, which treats the two-phase 
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flow of steam and water as a uniform mixture. However, it 
usually simulates the flow in the secondary loop with area- 
averaged variables. Given that coolant temperature varies 
significantly from the top to bottom of the U-tube bundles, 
both temperature and heat transfer coefficients vary consid- 
erably in the two flows. To examine the thermal-hydraulic 
characteristics of NPP SGs thoroughly, a code with a detailed 
model shall be developed. 

In this paper, we present a thermal hydraulic code for NPP 
SGs in a geometric model composed of the primary and sec- 
ondary loops, U-tube, and steam room. The unique secondary 
loop model is divided into hot and cold sides, and the flow in 
it is simulated using the four-equation drift flux model and 
is analyzed thermal-hydraulically through coupling with heat 
transfer of the tube wall. Finally, the code is used to imple- 
ment and verify the Qinshan NPP SG. The thermal hydraulic 
parameters are computed at 100%, 75%, 50%, 30%, and 15% 
power rates. 


Il. GEOMETRIC MODEL 


Given the complicated actual structure, the geometry of the 
nuclear SG should be simplified in modeling and thermal hy- 
draulic simulation. This work considers the U-tube NPP SG. 
The primary loop of the SG is assumed as a straight tube of 
equal length. The secondary loop is circular and consists of 
the water supply chamber and the descending and ascending 
channels. The ascending channel consists of the sub-cooling, 
boiling, and ascent sections. In the secondary loop, the divi- 
sion of hot and cold sides is defined by the flow direction in 
the primary loop. The side with the primary inlet flow (the 
hot side) is hotter than the side with the primary outlet flow 
(the cold side). The SG structure is composed of the primary 
and secondary loops, heat transfer tube, and steam room, as 
shown in Fig. 1. The straight section of the two loops is 7-m 
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Fig. 1. Simplified geometric frame of the SG. 


j 


long. The U-tubes are 0.022 m in diameter, with a total length 
of 16 m. 


HI. FIELD EQUATIONS 


In the U-tube SG, the flow types are of the single- and 
two-phase regions. Specifically, the flow in the primary loop 
and in sub-cooling sections of the secondary loop remains 
single-phase, whereas the two-phase region is illustrated by 
the flow in the boiling sections of the secondary loop. The 
two-phase regions are complicated in terms of flow and heat 
transfer. Thus, two sets of governing equations are estab- 
lished to model the single- and two-phase flows. 


A. Balance equations for single-phase regions 


In the sub-cooling sections of the U-tube SG, the single- 
phase region covers the primary loop, descending channel, 
and sub-cooling sections of the secondary loop. These re- 
gions utilize the single-phase flow model, and the respective 
balance equations of mass, energy, and momentum are as fol- 
lows: 


Op _ Apu) _ 
z + az = 0, (1) 
o ony P y 2 u SU 


Ahin- 1 es s/c BAI 
Cn nax Vo TERA 


Nucl. Sci. Tech. 25, 050601 (2014) 
ð ð, 
By Pu) tge Pe ae (3) 


B. Balance equations for two-phase regions 


The two-phase flow is mainly observed in the boiling sec- 
tion of the secondary loop of the U-tube SG. The equation of 
the four-equation drift flux model governs this region. This 
equation considers the velocity slip u, at the two-phase inter- 
face and the variation in void fraction along the flow path. In 
our work, the 1D, area-averaged governing equations of the 
four-equation drift flux model used are expressed as [8]: 
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where the subscript “m” pertains to fluid mixture parameters; 
Pm is density; Gm is mass rate; Um is velocity; Am is enthalpy; 
and uy = Ug — ur is the relative velocity of the liquid and 
gaseous phases. 


C. Heat transfer model 


We model the heat transfer between the primary and sec- 
ondary loops in the U-tube SG through the heat conduction 
of the tube wall. The heat transfer of the U-tube wall can 
then be simulated through 1D conduction in the cylindrical 
geometry. The convection heat transfer rate between the wall 
and the coolant is the source term, and the heat conduction 
equation is given by 


OT 1/0 OT 
Pep at È (E) a (8) 


r 


D. Correlations 


To close the field equations discussed above, we must de- 
termine the correlations in the thermal property of the fluid 
and wall materials, as well as the criteria for the different 
flow structures, resistances, heat, and mass transfers. The 
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unknown variables that must be derived from correlations in- 
clude pm, Cpm, Ig, he, hg, Twf, Twos Uwe, Uwe, Une, Une, q, Pts 
and pg. 

To identify the thermal property of water and steam, we ap- 
ply the formulas provided by the industrial standard IAPWS- 
IF97 [9]. The property of the Incoloy-800 alloy is consid- 
ered for the tube wall. Moreover, we apply the model of 
Taitel and Dukler in the structural criteria for flow [10]. Ac- 
cording to their model, bubble flow transitions to slug flow 
when bubble speed up is greater than Taylor bubble speed 
Ub given a low flow rate in the tube with a small diame- 
ter tube (Gm <2000 kg/(m? s)).This transition occurs when 
the void fraction is greater than 0.5 at an increased flow rate 
(Gm >3000kg/(m? s)), as shown by 


up = 1.53[9(pr — pg) / 070” 


uw = 0.35[gD(pr — pg) / p°”. 


The transition from slug flow to annular flow can be de- 
termined through the superficial velocity and the Kutateladze 
(Ku) number of the flow [11]. The transition is observed in 
the flow in the channels with small diameters when gaseous 
superficial velocity exceeds the critical superficial velocity 
Je,crit- However, the transition is initiated when the gaseous 
Ku number is greater than the critical Ku number, as ex- 
pressed by 


? 


(9) 
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The flow resistance in the U-tube SG considers the resis- 
tance to both gravitation and friction. The Darcy formula 
is applied in relation to the friction resistance of the single- 
phase flow. The split-phase friction model of Martinelli is 
employed [11] in relation to the friction resistance in the two- 
phase flow as: 


(10) 
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The convective heat transfer coefficient of the single-phase 
flow is calculated using the D—B formula with regard to flow 
in the primary loop and in the pre-heating section of the sec- 
ondary loop. The D-B formula is calculated according to 
Chen’s equation for the boiling section of the secondary loop 
[12], 
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The onset of nucleate boiling is computed using the model 
developed by Bergles and Rohsenow [13]: 


don = 0.015p" 156 (1.8A Tan), (13) 
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IV. NUMERICAL SCHEMES 
A. Numerical schemes of the flow field 


In our solution, we apply the semi-implicit difference 
scheme. We treat the convection terms in the mass and en- 
ergy equations, the pressure gradient, and the two-phase mass 
transfer in the momentum equation implicitly, whereas all 
other differential terms are examined explicitly. The stag- 
gered grids are applied in discretization, and two groups of 
control volumes are established in the same flow channel. 
The control volumes for pressure, void, density, and enthalpy 
are arranged in a staggered formation along with those for ve- 
locity. The mass and energy equations are discretized given 
the control volume groups i—1, 7, and ¿+1, which are also im- 
plemented by control volume groups j — 1, j, and j + 1 given 
the momentum equation. The values of the flow parameter 
are presumably uniform in all control volumes. Fig. 2 depicts 
the established staggered grids and control volumes. 


CV of mass and energy . 


a, h, p 
n, 


CV of momentum 


Fig. 2. Staggered grids for the discretization of flow conservation 
equations. 


In relation to the four-equation drift flux model used in the 
secondary loop, the semi-implicit discretization equations are 
listed below [13]: 
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The discretization equations of the single-phase flow in 
the primary loop are similar in form to those given above. 
To solve these discretization equations, we adopt a velocity— 
pressure correction scheme. First, the unknown pressure of 
the new time step is assigned a value equal to that of the old 
time step. Subsequently, the momentum equation is solved to 
estimate the velocity value of the new time step te Once 
the mixture mass, gaseous mass, and mixture energy equa- 
tions are rearranged, we obtain the following matrix equations 
for h?+1, att! and p?*?. 


m,t > 4 ? a 


n+1 n 
Au Ajo Ais Ani — hmi 


A21 A22 A23 arth gm | = 

Ası Asz Ass | | ppt! — pp (18) 
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By solving Eq. (18), we determine the pressure at the new 
time step ET. The new pressure is then used to correct ten- 
tative velocity and to obtain the final velocity at the new time 


n 
uag = at + ($2) (i-p). a9) 


J 


We apply a large time step, such as At = 10° s, for the 
steady state analysis. In this case, the time-derivative term is 
very small and can be disregarded in the discretization equa- 
tions. Thus, the balance equations above can then be applied 
to the steady-state solution. 


B. Solution for the U-tube wall conduction 


Equation (8) is integrated into cylindrical volume 27trdrdl 
at the time step At in relation to the heat conduction of the 
U-tube wall to generate the discretization equation for wall 
temperature. 
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C. Cross-iteration of heat balance and the natural circulation 
condition 


In the U-tube SG, heat transfer is simultaneous in the pri- 
mary and secondary loops; thus, flow and heat transfer in both 
loops must be coupled for solving. We adopt a coupled it- 
eration scheme that converges these factors when both heat 
balance and the natural circulation condition are satisfied. In 
heat balance, the heat transfer in the primary loop is equal to 
that in the secondary loop. In the natural circulation condi- 
tion, the head of driving pressure must meet the total pressure 
drop of the entire system, that is, 


DH=DP=)_AP;+5)_(AP,+5 AP. (21) 


In the iteration of heat balance, the heat fluxes in the pri- 
mary and secondary loops are initially assumed to be a group 
of values. Subsequently, matrix Eq. (18) is solved to deter- 
mine the flow parameters. The heat fluxes in the two loops 
are then computed in turn. In addition, the heat balance con- 
dition is validated. If the difference in heat flux between the 
primary and secondary loops is greater than a preset limit, 
the temperature of the coolant that enters the primary loop is 
corrected and a new iteration of heat balance is initiated. 

In the iteration of natural circulation, the dichotomy 
scheme is applied. First, the value of flow rate W is set, and 
the difference in pressure head and resistance is computed as 
f(W) = DH — D. The flow rate is then modified slightly 
to flip the sign of f(W’). The value of the flow rate is up- 
dated by W"t! = (W + W’)/2. The corrective iteration of 
W continues until f(W) meets a pre-set limit. Based on the 
theoretical model above, we therefore develop a code for the 
steady-state thermal hydraulic simulation of the nuclear U- 
tube SG. A numerical scheme is also established using MAT- 
LAB software. 


V. SIMULATION RESULTS 


The steady-state thermal hydraulic characteristics of the 
SG in the Qinshan 300 MW PWR are investigated with re- 
spect to the thermal hydraulic code presented for nuclear 
U-tube SGs. The grid gap measures 1.2m along the U- 
tube length. Moreover, this study considers five cases un- 
der different power conditions, namely, 100%, 75%, 50%, 
30%, and 15%. The results at the 100% power condition are 
compared with those simulated using the RELAP5 code [14]. 
Table | lists the required computation parameters given this 
power condition. 

Figures 3—7 show the computed steady-state thermal hy- 
draulic parameters at the 100% power level. The tube lengths 
are 0—8 m and 8-16m for the hot and cold sides of the sec- 
ondary loop, respectively. The results of the primary loop 
are plotted according to full tube length, whereas those of the 
secondary loop are plotted based on half tube length. Fig. 3 
presents the temperatures of the coolant in the primary and 
secondary loops and of the tube wall. The coolant tempera- 
ture decreases along the tube in the primary loop; in the sec- 
ondary loop, however, the inlet coolant is slightly sub-cooled. 
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TABLE 1. Condition parameters for 100% power. 


Parameters Value 
Flow rate of the primary loop (kg/s) 3333.3 
Inlet temperature of the primary loop (°C) 315.2 


Primary loop pressure (MPa) 15.3 
Flow rate of the secondary loop (kg/s) 259.86 
Inlet temperature of the secondary loop (°C) 215.6 
Steam room pressure (MPa) 5.43 
Water level of the steam room (m) 10.04 
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Fig. 3. Temperatures of the primary fluid (W, L1), secondary fluid (e, 
o),and U-tube wall(a, A) in the SG. The solid symbols represents the 
results of the current work, and the blank symbols denote the results 
obtained from RELAPS. 


Thus, the coolant temperature increases to saturation level af- 
ter a short distance. 

The simulation results with our code differ only slightly 
from those obtained with RELAPS. With respect to the 
coolant temperature of the secondary loop, our results are 
slightly lower than those derived from RELAP5. This may 
be attributed to different correlations assigned to the con- 
vection coefficient. In RELAP5, a modified correlation 
of the convection heat transfer coefficient (Nu = 2.0 + 
0.74Re'/? Prt/3) is applied to the single-phase liquid and 
sub-cooled boiling regions [14], whereas our study utilizes 
the D-B correlation. So, our technique generates a convec- 
tion coefficient value that is smaller than that obtained with 
RELAPS5. Temperature of the U-tube wall varies along the 
lengths in a manner that is almost similar to the coolant in 
the primary loop. This is ascribed to the fact that the heat re- 
sistance of the primary loop is much smaller than that of the 
secondary loop because the latter displays a noticeable foul- 
ing resistance. 

Figure 4 displays the gaseous and liquid velocities of the 
flow in the secondary loop under the 100% power rate. Both 
velocities increase continually in the secondary loop from the 
lower room to the steam room with tube heating and coolant 


Gas velocity (m/s) 


o 2 4 6 8 
Tube length (m) 


Liquid velocity (m/s) 


0 2 4 6 8 
Tube length (m) 


Fig. 4. Phase velocity in the hot (Mf, C) and cold((a, A) channels 
of the secondary loop under the 100% power condition. The solid 
symbols represent the results of this work, and the blank symbols 
denote the results obtained from RELAPS. 
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Fig. 5. Enthalpy of the fluid in the primary and secondary loops. 


boiling. The two-phase velocities increase in the steam room 
as a result of the expanding area. Furthermore, gas velocity 
is always higher than that of liquid because gas phase flow is 
affected by buoyancy. Nonetheless, the RELAPS results are 
5% higher than those of our code. 

Figure 5 exhibits the variation in coolant enthalpy along 
the tube lengths in the primary and secondary loops. Fig. 6 
shows the heat flux on the interior of the U-tube wall, which 
is equal to that of the exterior of the U-tube in steady-state 
analysis. The coolant enthalpy in the primary loop continues 
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Fig. 6. Heat flux on the interior of the U-tube wall. 
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Fig. 7. (a) Void fraction of the secondary fluid and (b) heat transfer 
coefficient at the tube wall under the100% power condition. 


to decrease along the tube length with heat transfer from the 
primary to the secondary loops, whereas that in the secondary 
loop continually increases throughout the process as depicted 
in Fig. 5. The heat flux on the interior of the U-tube decreases 
with tube length as the temperature difference between the 
tube wall and the coolant decreases along the tube (Fig. 6). 
Figure 7 displays the void fraction and the heat transfer 
coefficient along the tube lengths in the secondary loop. The 
coolant void fraction is higher in the hot side of secondary 
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Fig. 8. Fluid pressure in the primary and secondary loops under the 
100% power. 
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Fig. 9. Fluid temperature in the primary and secondary loops under 
the 100% power rate. 


loop than that in the cold side given that the heat flux in the 
hot side is higher. However, the coolant void fractions that 
enter the steam room from both sides of the secondary loop 
are similar as a result of lateral mixing. 

Figure 8 depicts the variation in the pressure of the primary 
and secondary loops. The pressure of the primary loop con- 
tinually decreases in the ascending part but increases in the 
descending part with the increase in gravitational potential 
energy. The pressures are similar at both sides of the sec- 
ondary loop and continue to progress downward along the 
tube length. 

Figure 9 displays coolant temperatures at five power rates 
in the primary and secondary loops. In the primary loop, in- 
let, outlet, and average coolant temperatures increase with the 
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Fig. 10. Variation in the circulation ratio and in circulation flow with 
power rate. 
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Fig. 11. Void fractions inthe cold and hot channels of the secondary 
loop given different power. 


increase in power rate. The temperature of the inlet coolant 
increases more quickly than that of the outlet coolant. Hence, 
the variation amplitude of temperature in the primary loop 
increases with power rate. Fig. 10 indicates that the satura- 
tion temperature of the coolant decreases when power rate 
increases. This finding suggests that the cooling capability of 
the secondary loop has been strengthened. 


Figure 10 presents the variations in circulation ratio and in 
circulation flow with power rate with regard to the SG. With 
the increase in power rate from 15% to 100%, the SG cir- 
culation flow initially increases at the small power rate but 
decreases when the power rate exceeds 50%. This result is 
induced by the coupled effect of driving pressure and circu- 
lation resistance. As the boiling length in the secondary loop 
increases with increasing power rate, the void fraction and 
driving pressure increase as well. Moreover, the circulation 
resistance increases with increasing flow rate; hence, the cir- 
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Fig. 12. Enthalpy in the (a) primary and (b) secondary loops under 
different power rates. 


culation flows downward along the tube length. The circu- 
lation ratio continually increases with power rate, as shown 
in Fig. 9. Furthermore, the mass flow of the vapor in the SG 
continually increases. 


Figure 11 shows the gaseous void fractions in both sides of 
the secondary loop at 75%, 50%, 30% and 15% power rates. 
This fraction increases along the tube lengths of both sides of 
the secondary loop. In addition, the gaseous void fraction is 
higher in the hot side than that in the cold side because boiling 
length is longer in the former. 


Figure 12 depicts the variation in coolant enthalpy with 
tube length in the primary and secondary loops under15%- 
75% power rate. The coolant temperature in primary loop de- 
creases along the tube length, where as that in the secondary 
loop is maximized. The enthalpy variation between the inlet 
and outlet of the two loops increases with high power rate. 


050601-7 


ZHANG Xiao- Ying et al. 


This result proves that the heat transfer process from the pri- 
mary loop to the secondary loop is strengthened. In the sec- 
ondary loop, coolant enthalpy increases more in the hot side 
than in the cold side. 


VI. CONCLUSION 


This study presents a steady-state thermal hydraulic code 
that was developed to thoroughly investigate the thermal hy- 
draulic characteristics of the nuclear U-tube SG. This code is 
based on the two-zone geometry model of secondary loop. 
Thermal hydraulic analysis was conducted using the four- 
equation flux model, and a cross-iteration solution was es- 
tablished to meet the conditions of heat balance and natural 
circulation. This solution is based on the staggered grids and 
the first-order scheme of explicit—implicit difference. The 
steady state thermal hydraulic characteristics of the SG were 
thus identified using the developed code for the QINSHAN 
I PWR under 100%, 75%, 50%, 30% and 15% power rates. 
Moreover, some important thermal and hydraulic parameters 
were identified for the primary and secondary loops. The re- 
sults obtained under the 100% power rate agree well with the 
results simulated using RELAPS. Hence, the established the- 
oretical model and numerical scheme can guide the design 
and safe operation of a nuclear U-tube SG. 


SYMBOL LIST 


p density, kg/m?; 

u velocity, m/s; 

t time, s; 

z distance, m; 

h enthalpy, kJ/(kg K); 

H heat transfer coefficient, W/(m? K); 
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